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ABSTRACT 

We present a family of spherical models for elliptical galaxies and bulges consisting of a 
stellar component and a central black hole. All models in this family share the same stellar 
density profile, which has a steep central cusp. The gravitational potential of each models is a 
linear combination of the potential generated selfconsistently by the stars and the potential of 
a central black hole. The relative importance of these two contributions is a free parameter in 
the models. Assuming an isotropic dynamical structure, almost all kinematical properties of 
these models can be calculated analytically. In particular, they form the first simple dynamical 
models for galaxies with a central black hole where the distribution function and differen- 
tial energy distribution can be written completely in terms of elementary functions only. We 
also present various extensions of this family to models with anisotropic orbital structures. 
Also for these dynamical models, the distribution function and its moments can be expressed 
completely in terms of elementary functions. 

This family is useful for a large range of applications, in particular to generate initial 
conditions for N-body and hydrodynamical simulations to model galactic nuclei with a central 
black hole. 

Key words: black hole physics - celestial mechanics, stellar dynamics - galaxies: kinematics 
and dynamics - galaxies; structure 



1 INTRODUCTION 

During the past few years, various numerical dynamical modelling techniques have been developed to construct accurate dynamical models 
for galaxies (Dejonghe 1989; Emsellem, Monet & Bacon 1994; Rix et al. 1997; van der Marel et al. 1998; Gerhard et al. 1998; Cretton et 
al. 1999; Gebhardt et al. 2000a; Verolme & de Zeeuw 2002). The state-of-the-art techniques can handle various degrees of complexity such as 
the capability to include non-trivial geometries and higher-order kinematical information. Nevertheless, analytical models remain interesting 
and important for various reasons. They can provide a simple testbed where various physical processes can be investigated, or where new 
modelling or data analysis techniques can be explored. Analytical models are particularly useful to generate the initial conditions for N-body, 
SPH or Monte Carlo simulations. Most attention has been devoted to the construction of spherical self consistent models, i.e. models in which 
the stars move in a potential generated completely by the stars themselves. Famous examples include the Plummer model (Plummer 1911; 
Dejonghe 1987), the isochrone sphere (Henon 1959, 1960) and the Jaffe sphere (Jaffe 1983) and the Hernquist model (Hemquist 1990; Baes 
& Dejonghe 2002). Such models often serve as a template model for elliptical galaxies, globular clusters or the bulges of spiral galaxies. 

Recent observational evidence indicates, however, that the central regions of these stellar systems cannot always be faithfully represented 
by selfconsistent models. The existence of supermassive black holes in the nuclei of galaxies has been suspected for a long time, as accretion 
onto massive compact objects was regarded as the only reasonable explanation for the existence of active galaxies and quasars (Salpeter 1964; 
Soltan 1982). During the last decade, quiescent supermassive black holes have been detected in the centre of the Milky Way (Ghez et al. 2000; 
Schodel et al. 2002) and in the nuclei of virtually all external galaxies which are nearby enough to resolve the sphere of influence of the black 
hole (e.g. see list in Tremaine et al. 2002). The masses of these black holes are roughly between a million and a few billion solar masses 
and are tightly coupled to structural parameters of the host galaxies (Kormendy & Richstone 1995; Ferrarese & Merritt 2000; Gebhardt et 
al. 2000b; Graham et al. 2001; McLure & Dunlop 2002; Ferrarese 2002; Baes et al. 2003; Marconi & Hunt 2003). Recently, evidence for 
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the presence of intermediate mass black holes in the centre of globular clusters has been reported (Gebhardt, Rich & Ho 2002; Gerssen 
et al. 2002, 2003), although this evidence is still controversial (Baumgardt et al. 2003a,b; McNamara, Harrison & Anderson 2003). These 
observations clearly indicate that there is a need for simple but realistic analytical dynamical models with a central black hole. 

The first requirement for the construction of models with a black hole is a potential-density pair with a cuspy central density profile. 
Tremaine et al. (1994) demonstrate that spherical isotropic systems need a central density profile that diverges at least as r"^''^ to be able 
to support a central black hole. The need for such cuspy models also follows from observations: HST photometry has shown that ellipticals 
and the bulges of spiral galaxies generally have central density cusps that diverge as at small radii with 0^7^2.5 (Lauer et al. 1995; 
Gebhardt et al. 1996; Ravindranath et al. 2002; Seigar et al. 2002). Also for the Milky Way there is evidence for a central density cusp in 
the innermost regions (Genzel et al. 2003). A very useful one-parameter family of spherical models with this characteristic, which we refer 
to as the 7-models, was introduced independently by Dehnen (1993) and Tremaine et al. (1994). This family has a simple analytical density 
profile which diverges as r~'^ in the central regions (0 < 7 < 3), and includes the Hemquist and Jaffe models as special cases. Zhao (1996) 
generalized this family further to a very general three-parameter family of models, the so-called (a, (3, 7)-models. For many of these models, 
interesting dynamical properties such as the intrinsic and projected velocity dispersions, the distribution function and the differential energy 
distribution can be calculated analytically if one assumes an isotropic selfconsistent dynamical structure. 

The second step in the construction of dynamical models with a cuspy core is to add to the potential of the selfconsistent model an extra 
contribution from the black hole, and re-calculate the dynamical properties with this new potential. The calculation of the (intrinsic and/or 
projected) velocity dispersions in the presence of a black hole is not very difficult, and can usually be performed analytically for those models 
where the dispersions can be calculated analytically in the selfconsistent case. For example, for the sets of potential-density pairs considered 
by Tremaine et al. (1994) and Zhao (1996), the addition of a black hole was not a problem for the calculation of the velocity dispersion profile. 
The reason is that the intrinsic and projected velocity dispersions are just Unear functions of the potential, and therefore linear functions of the 
black hole mass. Many other interesting kinematical properties, in particular the phase space distribution function, depend on the potential 
in a strongly non-linear way, however. The construction of dynamical models in which these more complicated kinematical quantities can be 
expressed analytically in the presence of a black hole proves to be more difficult. Apart from the asymptotic behaviour, these characteristics 
usually need to be calculated numerically (Tremaine et al. 1994; Dejonghe et al. 1999). 

The work by Ciotti (1996) provides a first attempt to construct completely analytical dynamical models for galaxies with a central black 
hole. In an effort to construct realistic dynamical models for elliptical galaxies embedded in massive cuspy dark matter haloes, he considers 
a set of two-component Hernquist models. They consist of a stellar component and a halo component, both of which are modelled as a 
Hemquist profile, but with different masses and core radii. The interesting aspect of his work for our goal is that the masses and core radii of 
each component can be taken arbitrarily. Setting the core radius of the halo component to zero, his two-component model degenerates into a 
Hernquist model with a central black hole. Ciotti demonstrates that the distribution function and differential energy distribution of Hernquist 
models with a central black hole can be expressed analytically, albeit as rather cumbersome expressions. For example, the distribution function 
can be written as the derivative of a compUcated function involving incomplete elliptic integrals and Jacobian functions. Nevertheless, his 
work presents the first model for galaxies with a central black hole where most of the kinematical properties can be calculated analytically. 

In this paper, we present a detailed kinematical analysis of a family of spherical models for elliptical galaxies and bulges consisting of 
a stellar component and a central black hole. In Section 2 we define our family of models. AU models in this one-parameter family share the 
same stellar density profile, which corresponds to one of 7-models introduced by Dehnen (1993) and Tremaine et al. (1994). The potential 
of the models is a linear combination of the potential generated by the stars and the potential of a central black hole. The importance of 
both components to the total potential can be set arbitrarily by varying the parameter n, which represents the mass of the central black hole 
relative to the total mass. The reason why we chose this particular family of models is that most of the interesting kinematical properties 
can be expressed completely analytically for all values of jx. In Section 3 we describe some intrinsic kinematical properties of this family of 
models, such as the distribution function, the differential energy distribution and the moments of the distribution function. In particular, we 
give closed analytical expressions for these quantities, completely in terms of elementary functions. In Section 4 we describe some observed 
kinematical properties as they are projected on the plane of the sky. Most of these quantities can be expressed analytically involving simple 
incomplete elliptic integrals. In Section 5 we present a number of generalizations to models with a anisotropic orbital structure. In particular, 
we discuss models with a constant anisotropy profile, models with a distribution function of the Osipkov-Merritt type and models with a 
completely radial orbital structure. Also for these models, the distribution function and its moments can be calculated completely analytically 
for all values of the central black hole mass. Finally, a discussion is given in Section 5. 



2 DEFINITION OF THE MODELS 



The starting point for our family of models is the stellar luminosity density profile 

L* fr\ -5/2 / r 



Li, rr\ -5/2 / r\ -3/2 

Pir) = ^{A (1 + -) , (1) 
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where L* is the total stellar luminosity of the system and c is a scale length. The (positive) potential generated by this stellar distribution is 
easily found through Poisson's equation, 



T / N ^ r ^ 2GAf, r + c ^ . 

*-M = -**M = ^— I y^-ij , (2) 

where is the total stellar mass. The isotropic selfconsistent model with this potential-density pair is a particular case from the family of 
7-models considered by Dehnen (1993) and Tremaine et al. (1994), corresponding to 7 = | or 7; = i. A number of kinematical properties 
for this selfconsistent model, such as the velocity dispersion and distribution function has been derived in these papers. Our aim is to consider 
systems consisting of two components: a stellar component and a central supermassive black hole. In order to find the total potential of such 
a two-component system, we need to add to the stellar potential an extra contribution from the black hole. We model the central black hole 
as a point potential with mass M,, such that we obtain 



*W = ^ (3) 



In the remainder of this paper, we will work in dimensionless units; we set the gravitational constant G, the scalelength c, the luminosity I/* 
and the total mass M^, + M, of the model equal to unity. We introduce the parameter /i as the fraction of the black hole mass to the total 
mass of the galaxy. With these conventions, we can write the potential-density pair of our model as 

^('■) = ^ r5/2(l + r)3/2' 



(5) 



In this paper, we consider isotropic kinematical models implied by the potential-density pair formed by equations J4} and Js}. This potential- 
density pair defines a one-parameter family of isotropic dynamical models, where is a free parameter that can assume any value between 
and I. The models corresponding to these values of /i have a particular meaning: the former corresponds to a selfconsistent model without a 
central black hole, whereas the latter represents a galaxy with a central black hole and a negligible stellar mass-to-light ratio. 

As an important remark, the reader should note that the convention we use is different from the convention used in e.g. Tremaine et 
al. (1994) and Zhao (1996). In these papers, /i denotes the black hole mass relative to the stellar mass, and the normalization is such that the 
stellar mass is set to unity. We prefer to set the total mass of the galaxy equal to unity however, because (1) all models then have the same 
behaviour at large radii and (2) this allows us to study the entire range of models from selfconsistent models without black hole to models 
completely dominated by the black hole potential. 



3 INTRINSIC PROPERTIES OF THE MODELS 
3.1 Basic properties 

The total mass density, circular velocity curve and the cumulative mass function are readily obtained 
1 - /i 1 



ptot(r) = 



2, . 1 - M 

(r) = 



8tt r5/2(i + r)3/2 



+ 



— + ^^. 

+ r 

The isotropic velocity dispersion a — a,. — as = o^ can be obtained from the solution of the Jeans equation. 



M{r) = (1 - m) 



M{r) p{r)dr 



(6) 
(7) 

(8) 
(9) 



p{r) 

The cumulative mass function, which depends linearly on the black hole mass, is the only /^-dependent quantity in this expression. Con- 
sequently we can write the velocity dispersion as the sum of separate contributions from the stellar mass and the black hole. The result 
is 

jl + r 



a (r) 



1 - 2r + 6r + 12r'' - 12r^(l +r)ln 



1 



+ 



2p. f 1 + r 
35 



5 - 8r + 16r^ - 64r' ~ 128r* + USr'^^^VTTr 



At large radii, the velocity dispersion goes to zero as 



12r 



+ ■ 



(10) 



(11) 
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Figure 1. Some intrinsic and projected kinematical properties of tlie model for various values of the parameter /i. The first panel shows the intrinsic velocity 
dispersion profile u{r) [solid lines] and the projected velocity dispersion profile crp(R) [dotted lines], the second panel shows the distribution function f{£) 
and in the right panel the differential energy distribution M{£) is plotted. The colour code is the same for all panels: = (black), ^ = 0.01 (yellow), 
fi = 0.1 (red), /i = 0.5 (green) and fJ, = 1 (magenta). 



At small radii, the velocity dispersion diverges as 

In the left panel of figureQwe plot the intrinsic velocity dispersion profile for various black hole masses. Without a central black hole, the 
velocity dispersion (and the circular velocity) have a weak r^^* divergence at small radii. When a black hole is present, these divergences 
become stronger and both quantities diverge as r^^^'^ . 



3.2 The distribution function 



The ultimate goal of the dynamical modeller is to obtain an expression for the phase space distribution function /(r, v), which contains all 
kinematical information on a system. For non-rotating isotropic spherically symmetric systems, the phase space distribution function depends 
on the position and velocity vectors only through the binding energy £ — ii){r) — \v^, which is a integral of motion. The distribution function 
/(r , v) — f{£) is thus a one-dimensional rather than a six-dimensional function. The key to calculate the distribution function corresponding 
to a given potential-density pair is the augmented density p{^), i.e. the density written in terms of the potential. With this function, we can 
calculate the distribution function through the Eddington (1916) relation, 



fin = 



1 



d?p d* 



1 

7f 



dp 

d* 



(13) 



,0 d*2 VC V"-/*=0J 

The second term in this expression vanishes for the family of models we consider in this paper, both with and without black hole, because 
oc at large radii. To find an expression for the augmented density, we need to invert the potential to a relation r(^) and determine 
p(*)=p(r(*)). 

We first consider the special case when the galaxy contains no central black hole (p — 0). The augmented density is readily found, 
1 *4 (4 + *)4 



2567r (2 + *)3 ' 
and the resulting distribution function is (Tremaine et al. 1994) 



(14) 



fin 



V2 
8967r3 



42 (3 - 32£ - 8£^) . , IT (63 - 693£ - SGTOf^ - 7410£:^ - 4488£:'' - 1448£:^ - 240£:'' - 16£'^)V£ 

; ::;T7r7:; arcsinh \ -. -—. 

(2 + £:)9/2 V 2 {2 + £)'>' 



(15) 



The behaviour at small binding energies is 

OTV^ \ 7 

whereas at large binding energies we obtain 



1-^ + - 



fin^^ 



%/2 „7/2 



567r3 



(16) 



(17) 



Another special case is the other side of the spectrum (/i = 1), which corresponds to a galaxy where the black hole completely dominates 
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the potential. In this case the augmented density and distribution function read 
1 

" 8^ (l + *)3/2' 



fin 



V2 
1287r3 



(l+f)3 



3 (3 — 5iS) arctan VS 



The behaviour at small binding energies is 



2^2 



^5/2 



1 - 



15£ 
7 



whereas at large binding energies we obtain 



/(£) 



15^2 



£ 



(18) 
(19) 

(20) 

(21) 



2567r2 

For the general case < ^ < 1, i.e. when both the stars and the central black hole contribute to the potential of the galaxy, the calculation of 
the distribution function is more tedious. In order to find a convenient expression for the augmented density, we write expression J4j as 



on 



with 



a = a(r) = 



1 



A relation a (5') can be found by inverting the expression 
^{a) = 2 (1 - At) (a - 1) + M (a^ - 1) , 
which yields 



a(*) = 



Ai - 1 + VI + 



(22) 



(23) 



(24) 



To calculate a closed expression for the distribution function, we can now take different approaches. The first approach is a direct application 
of the Eddington formula <13t . When we substitute the explicit expression for into the formula 122\ . we find after some algebra that the 
second derivative of the augmented density can be written in the form 



-(*) - 



327rAi3 (2-Ai + *)^ 



J'2(*,m) 



l + /iV[')3/2 



(25) 



where Pi and P2 are polynomials with integer coefficients in 9 and jj.. This form is suitable for analytical integration in the Eddington 
relation: the two parts of the integrand are basically the combination of a rational function of 4' and the square root of first and second order 
polynomials in \&. Such a function can be integrated analytically and the resulting integral can be written in terms of elementary functions 
only. The distribution function can also be calculated in a more elegant way by transforming the Eddington relation <13t to an integration 
with respect to a. We obtain 

da 



fin-- 

with 

as = 



1 





■dp 


/ d*' 


Ji da 


da/ 


da 



V'£-2(1-At)(a-1)-M(a2-1) 



1 + yi + m£ 



With the expressions <22t and <23> the first factor in the integrand becomes 



_d_ 

da 



dp 
da / 



d* 
da 



1 (a^ 
16tv 



if [12(1 - m) + ^5^^a + 16(1 - p)a^ + IS/xa^ + 20(1 - /x)a'^ + IS/xa^ 
{1 — fj, + /ia)2 



(26) 



(27) 



(28) 



This function can be decomposed in partial fractions, and each of the resulting integrals can be calculated analytically in terms of elementary 
functions. Using either of both approaches, we obtain that the resulting distribution function can be written explicitly as 



fin = 



V2 
1287r3 



Ai3 (2- + £)*{! + ^i£) 



+ 



+ 



3 (25 - 40ai + I2fi'^ + 5fj.£) 



,7/2 



arctan ^/ p£ 
3 (1 - At) (3 - S/i + 12aj^ - 32g - Slfi£ - 8£^) 



(2-p + f)9/2 



arctanh 



^{2^fi + £)£ 
l + £ 



(29) 



where X{£,ii) is a polynomial with integer coefficients 
X{£,ii) = 3(400 - 1440/1 + 2072/1^ 



1541/1^ + 622/1* 



116/1^) + (2400 - 6400/1 - 5352/i^ - 503/i^ - 1658/i* + 1002/i^ + 220/1^) £ 

2175/ + 373/1^) + 2 (300 + 150/1 - 1356^1^ + 1116/i^ - 253p'') f ^ 

16/l)£^ (30) 

The first two terms in the distribution function <29> both diverge as /i^^ when /i approaches zero, but the diverging terms of course cancel 
out, such that the distribution function reduces to the expression <15> in the limit p ^ 0, as required. On the other hand, it is straightforward 



+ (1800 - 2600/t - 1486/1^ + 4030/i^ 
+ (75 + 400p - 864/t^ + 328p^) £^ + 5p (13 
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to check that the distribution function \29\ reduces to the expression <19t in the limit ^ 
function <29> behaves as 

Stt'^ [ 7 
and at large binding energies we have the expansion 



1. At small binding energies, the distribution 



(31) 



15\/2 £ 



32(1 - ^) ^ 25 - 40^ + 12/i2 



+ 



5^8 



+ ■ 



(32) 



In the second panel of figure Q we plot the distribution function for different values of the black hole mass. As the central potential well 
is infinitely deep, orbits with all binding energies are allowed. Both without and with a central black hole, the distribution function is 
a monotonically rising function, which guarantees the stability of the model against radial and non-radial perturbations. At small binding 
energies the influence of a black hole is negligible. This is also visible in the asymptotic expansion <31> . where the leading term is independent 
of the black hole mass fraction fi. The same leading terms is then obviously found in the asymptotic expansions <16> and <20> . At large binding 
energies, i.e. in the central regions of the system, the influence of the black hole is important. The distribution function becomes less steep 
in the large binding energy limit when a black hole is present: its slope changes from | to 1. The larger the black hole mass, the smaller the 
values of the distribution function. 



3.3 The differential energy distribution 



The distribution function f{r, v) = f{£) contains all available kinematical information on the system, but is rather hard to interpret. In par- 
ticular, f(£) does not represent the number of stars per unit binding energy. The quantity that describes this useful characteristic is the differ- 
ential energy distribution A/'(f). For an isotropic spherical system, the differential energy distribution can be written as N[£) = f{£) g{£), 
where g{£) is called the density of states and represents the phase space volume per unit binding energy (Binney & Tremaine 1987). This 
function can be calculated through the equation 
dr 



g{£) = le^TT^ 



d* 



d*. 



(33) 



(34) 



For the selfconsistent model without black hole, the density of states can immediately be calculated 

" 1 _ 21 + 9£: + £'^ 

and the differential energy distribution is found by multiplying this expression with the distribution function <15> . At small binding energies 
the differential energy distribution asymptotically behaves as 



g{£) = V2^' 



N{£): 



12£ 



and at large binding energies we have the expansion 
3 



J^{£): 



2£^ 



(35) 



(36) 



If the potential of the galaxy is completely dominated by the black hole, we obtain the well-known and very simple expression for the density 
of states function, 

V2 7V'' 



g{£) 



£5/2 



(37) 



Combining this expression with equation <19> for the distribution function, we find that the differential energy distribution asymptotically 
behaves as 



in the small binding energy limit, whereas at large binding energies we have the asymptotic expansion 



U{£): 



IStt 
128£'3/2 



55 



+ ■ 



(38) 



(39) 



In the general case < ^ < 1, the calculation of the density of states is less straightforward. Similarly as for the distribution function, we 
have calculated the function g{£) in two ways. The first approach uses a brute force calculation of the definition <33> . A more subtle approach 
consists of transforming the integral in this expression to an integration with respect to a. We obtain 

a 



g{£) = 32%/27r^ 



(a2 _ 1) 



- \/2(l - At)(a - 1) + - 1) - da, 



(40) 
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where a£ as in equation Ml\ . This integral can be calculated analytically and results in 

2(1 - ^J) [24(1 - m) - 4(3 - 10^ + 4At^)£: - 2(9 - 11^)^^ - Sf^] ^ 



g{£) = %/2^' 



+ 



3£:2(4-4^i + £:)2 

1 - (1 - ^)£: + (1 - ^)£:^ 

(21 - 70^ + 80At^ 



I^TT — arctan \/ /xf 



32^^) + (9 - 19^ + 10^^)£: + (1 - 



arctan 



2/i - 1 



(41) 



(4-4/i + f)5/2 

In the limits /i ^ and /i ^ 1, this function reduces to the expressions <34t and <37t respectively. The differential energy distribution A/" (£) 
of the model we consider can thus be written completely in terms of elementary functions for all values of ^. At small binding energies, the 
differential energy distribution behaves as 



3(4 + ij)£ 
7 



and at large binding energies we find the expansion 



(42) 



100-215^+112/1^ 8192(1--^)^ 



(43) 



^ ' 128£:»/2 [ 57r^ VS \ 5m 457r2 /i 

In the right-hand panel of figure Q we plot the differential energy distribution for models with various black hole masses. Without a black 
hole, the differential energy distribution is a monotonically decreasing function that has a finite value at £ = and that decreases smoothly 
to zero as £~'^ at large binding energies. In spite of the strong divergence of the distribution function at large binding energies, the system 
thus has no stars at rest in the centre of the galaxy. Both of these properties (a finite value at f = and convergence to zero at large £) 
are common properties for all selfconsistent models in the family of 7-models (Dehnen 1993). Adding a supermassive black hole does not 
change the behaviour of the differential energy distribution very drastically. The differential energy distribution still reaches the same finite 
value I at f = 0, and then smoothly decreases for increasing binding energies. The slope of the differential energy distribution is a weakly 
decreasing function of /i. At large binding energies, the differential energy distribution still converges to zero, but the slope of the convergence 



is weakened from —2 to 



For increasing black hole mass, the differential energy distribution assumes smaller values at small binding 



energies and converges less rapidly to zero at large binding energies. Models with increasing black hole masses hence contain stars which 
are in the mean increasingly strongly bound to the galaxy. 



3.4 The moments of the distribution function 

It is useful to consider the moments of the distribution function with respect to the velocities. In a general non-rotating spherical system, the 
moments of the distribution function are defined as 

M2i,2j,2fc {r) = j I j f{r, v)vl'- vl^ vf dw, (44) 

where the integration covers the entire velocity space. For isotropic spherical systems the distribution function is symmetric in the velocities, 
and one defines the isotropic moments as 



/■V2^C'-) 

M2n(r) = 47r / f{r,v) 
Jo 



%/2*M 

dw. (45) 



It is easily shown that the relation between the general moments and the isotropic moments is 

I ^ 1 r(z + i)r(j + i)r(fc + i) 

M2.2„2.(r) = r(,^^^fc^3) M2(,+.+.)(r). (46) 

Knowledge of the isotropic moments is thus sufficient to determine all the moments. A direct integration of the definition is not the most 
obvious way to calculate the moments. A more interesting way is to use a formula that expresses the augmented moments, i.e. the moments 
written as a function of the potential, as a function of the augmented density (Dejonghe 1986), 

A2n(*) = £ (* - P (*') d*'. (47) 

For the model we consider, we can write this expression as 

~ f.TA ~ r\ 1 (2n + l)!! (l - a'")^ (1 - M + M"') ro n u '\ ^ I ^ '2^"-!^' r/is^ 
/i2n(4') = M2n(a) = — jj||-y ^ '—^^ [2 (1 - fi) [u - a ) + fi [a - a )\ da. (48) 

For all integer values of n, the integrand of this integral is a simple power series in a' and can easily be integrated analytically. The 2n'th 
moment can be written as a polynomial in /i of order n, where each coefficient is a sum of powers in a and terms in In a. The true isotropic 
moments fj,2n (r) can subsequently be found by substitution of a = ^y (1 + r)/r into the result. For example, for the second order moment, 
one obtains 



_ . . 3 
P2(a) = - 

TV 



(^-'^n-8^-l2-l'^«+ — -T + 24 +^ -a + 35-" + Y-y+28 



(49) 
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If we substitute a — + r) /r and /i2 = 'ipa^ into this equation, we recover the expression llOt . 



4 OBSERVED PROPERTIES OF THE MODEL 
4.1 The surface brightness 

The surface brightness profile can be found by projecting the luminosity density profile on the plane of the sky, 
I{R) = 2 



p{r) r dr 



Vr 



i?2 



(50) 



As all models in the family we investigate share the same stellar density profile, they will obviously share the same surface brightness profile 
It can be written in terms of incomplete elliptic integrals of the first and second kind, 

i + ^l-ii E^- 



I{R) = - 



(51) 



R+lJ ^R{R+1) R+1^ 

In the Appendix we give the relevant formulae to transform this formula (and similar formulae which will appear in the remainder of this 
discussion) to a more convenient form for R < 1 and R ~ 1. At large projected radii, the surface brightness decreases as 
1 4 



liR) 



+ ■ 



i6i?3 V '^R 

and the behaviour at small projected radii is 



liR) 



1 _L p-3/2 

2 r2 



3r^ 

167r2 



R 



(52) 



(53) 



where F = r(i) ~ 3.62561. The cumulative surface brightness profile can be found by integrating the surface density over a circular surface 
on the plane of the sky, 

S{R) = 2n I{R')R'dR'. (54) 
Jo 

A more convenient expression can be found by substitution of the definition <50> into this expression and partial integration, 

p oo 

S{R) = 1-4tt p{r)^r'^ - K^rdr. (55) 

Jr. 

For the luminosity density profile we consider, this expression can also be written in terms of incomplete elliptic integrals of the first and 
second kind 



S{R) = 2^/R{R+l)\ 



V^RTT I 4' 



The effective radius Rai is obtained by solving the equation S{Rai) = |; a numerical evaluation gives i?eff = 0.244955. 



(56) 



4.2 The projected velocity dispersion 

The projected velocity dispersion can be found by projecting the intrinsic velocity dispersion profile on the plane of the sky, 
p(r) (J^(r) r dr 



a'JR) 



I(R) 



or after substitution of the Jeans solution j9j and partial integration, 

a',{R) 2 



p(r)M{r) Vr^ - R^ dr 



(57) 



(58) 



I(R) JR 

Similarly as for the intrinsic dispersion, we can write the projected velocity dispersion as the sum of separate contributions from the stellar 
mass and the black hole. The contribution of the stellar mass leads to an integral which can be written in terms of elementary function, 
whereas the contribution from the black hole leads to a more complicated integral which can be written in terms of the incomplete elliptic 
integrals of the first and second kind, 

"4i?^-l 12i?^ + l i-m^ 



I{R)a;[R) = {l~^l) 



1057r 



5- 



127ri?2 



47r 



■X{R) 



'96R^ 2 - 17) VR + T / tt 
i?2 + R-i/2 ''^ I 4 ^ 



7?+ 1 



192i?'' - 82R^ - 5^1 TV 

:^=^= F — 

R^/^^/RTT \ 4' 



R+1 



(59) 



Dynamical models for galaxies and bulges with a black hole 9 



with X{R) a continuous real function defined as 



XiR) 



1 



arccosh 



: arccos — 



if < f, 



if i? > 1. 



y/W- 1 

At large projected radii, the projected velocity dispersion profile behaves as 
f57r^(4-/i) -512 



157ri? 



1287r J? 



whereas at small projected radii we obtain 



i_fi 

6V27r3/2 y/R 



167r2 



■R + 



847r2 R 



3(13447r'' - 5r* 



■R 



(60) 



(61) 



(62) 



In the first panel of figureQwe plot the projected velocity dispersion profile for various values of the black hole mass. Notice that the projected 
velocity dispersion profile is nearly similar to the intrinsic dispersion profile. At a given projected radius R, the projected dispersion is a 
weighted mean of the dispersion along the line of sight. As the density is a strongly decreasing function of radius, stars around the tangent 
point r ~ R will strongly dominate this mean, such that the projected velocity dispersion is nearly equal to the intrinsic dispersion at this 
point. As a result of this similarity, the effect of a black hole on the projected velocity dispersion is obviously similar as the effect on the 
intrinsic velocity dispersion. At large projected radii, there is obviously no effect from the black hole, whereas a black hole changes the slope 
of the projected dispersion profile at small projected radii from ~| to — i. 



5 EXTENSION TO MODELS WITH AN ANISOTROPIC ORBITAL STRUCTURE 

The models in the family we have described all have an isotropic dynamical structure, and consequently a distribution function f{£) that 
depends only on the binding energy. General anisotropic distribution function in spherically symmetric systems have distribution functions 
that depend on two integrals of motion, usually the binding energy and the modulus L of the angular momentum vector. The construc- 
tion of general anisotropic dynamical models f(£,L) that correspond to a given spherical potential-density pair is discussed in detail by 
Dejonghe (1986). Also for anisotropic models, the key ingredient for the calculation of the distribution function and its moments is an aug- 
mented density. For anisotropic models, the augmented density is an explicit function of both radius and potential, i.e. a function , r). For 
each potential-density pair, infinitely many of these augmented density functions and thus anisotropic dynamical models can be constructed. 
The only requirements are that the condition 

p{9{r),r) = p{r) (63) 

is satisfied and that the corresponding distribution function is non-negative over the entire phase space. Unfortunately the formulae to calculate 
the distribution function and its moments from the augmented density are significantly more complicated than the Eddington formula <13> 
in the isotropic case. A number of different approaches exist, including for example an approach with combined Laplace-Mellin transforms. 
There are a number of special cases, however, for which the construction of anisotropic dynamical models is not much more demanding 
than the construction of isotropic models. This has been illustrated by Baes & Dejonghe (2002), who constructed three different families 
of analytical dynamical models for the selfconsistent Hemquist potential-density pair, with widely different dynamical structure. In this 
section we shortly describe how our family of isotropic dynamical models with a central black hole can be generalized to models with a 
constant anisotropy, models with an Osipkov-Merritt type distribution function and models with a completely radial orbital structure. Further 
generalizations, such as Cuddeford (1991) models, are also possible. 



5.1 Models with a constant anisotropy 

A special family of dynamical models corresponds to models with a augmented density p{^, r) that is a power law of r, 

=pA(*)r-2'', (64) 

where /3 < 1 and where the function Pa(*I') is determined by the condition 163 > . It is well-known that augmented densities of the form J64t 
correspond to models with a constant anisotropy (Dejonghe 1986; Binney & Tremaine 1987). The distribution function is a power law of the 
angular momentum, and can be found through an Eddington-like formula 

f(F n- 2' f ^"P^ r6s^ 

' ' (27r)3/2 r(l-/3)r(i + /3) Jo d*2 (£:_^,) 1/2-/3 • 
To calculate the distribution function for the model defined by the potential-density pair 0-^5}, we can transform equation <65t to an 
integration with respect to a in the same way as in section l7!2l We directly obtain 

_ , , 1 (g^ - If-'^ ,,,, 
''^("^=8;^ • ^^^^ 
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If /3 is a half-integer number, the integrand in formula <65> is a rational function of a, and if (5 is an integer number, the integrand is the 
product of a rational function and a square root of a quadratic term in a. The distribution functions can consequently be written in terms of 
elementary functions only for all integer and half-integer values of (3. The radial and tangential velocity dispersions can be found through a 
formula similar to equation ^9} 

2. ^ <^iM '^li.r) 1 r°° M{r)pA{r)dr 

"^('^)= W-T3^ = ^X ~- ■ ^''^ 

These expressions can also be written in terms of elementary functions for all integer and half-integer values of f3. Similarly, analytical 
expressions can be derived for all higher-order moments if /3 is an integer or half-integer number. 



5.2 Osipkov-Merritt models 

Osipkov (1979) and Merritt (1985) developed an inversion technique for another special class of spherical anisotropic dynamical models, 
namely models where the distribution function depends on energy and angular momentum only through the combination Q = £ — /2r^, 
with Ta the so-called anisotropy radius and under the assumption that f{£,L) = for Q < 0. Such models are characterized by an augmented 
density of the form 

P(*,0=(l + J) (68) 

where the function Pq(\[') is determined by the condition J63> . The Osipkov-Merritt models have an anisotropy profile of the form /3(r) — 
/ (r^ +rf,), i.e. they are isotropic in the centre and become completely radially anisotropic at large radii. For the calculation of distribution 
functions for Osipkov-Merritt models, one can use a formula very similar to the Eddington formula. 

For the model we consider, we obtain 
1 fa^ - 1)2 

Pq(«) = p{a) + ^ ^ 3 ^ (70) 

The distribution function of the Osipkov-Merritt models can be calculated in a very similar way as the distribution function of the isotropic 
models. We find after some algebra 



f{£,L)^f{Q)' ^ 



1287r3 r| 



(2-m + 0)*(1 + mQ) 



3 (1 - (19 - 24At + 16p^ - leg + 23^Q - 4Q2) / ^{2- p + Q)Q 

3,rctciiiii f 



(2-At + Q)3/2 c......... ^ , , (71) 

where f(Q) represents the isotropic distribution function <29> and Xq{Q, p) is a polynomial with integer coefficients 

Aq(Q, p) = (199 - 383m + 264^2 - 80^*^ + 16/) + (67 + 58^ - imp^ + 96^i^) Q (14 + 7ai - bp^) + 2{l~ p) (72) 

The velocity dispersions for the Osipkov-Merritt models can be found through the formula 

2m = '^^'('^) = '^^(^) ^ _L_ r ^ir)pQ{r)dr 

PQWi, r2 • 

This integral can easily be evaluated analytically, and also the higher-order moments of the distribution function can be expressed completely 
in terms of elementary functions for all values of ra and p. 



5.3 Models with a completely radial orbital structure 

It deserves some special attention that the distribution function of both the models with constant anisotropy and the Osipkov-Merritt models 
corresponding to our potential-density pair l4j-|3 are positive for all values of /3 and ra and for all values of the relative black hole mass 
p. This is quite unusual. In the limit ra — > or /? ^ 1 for the Osipkov-Merritt and constant anisotropy models respectively, the entire 
galaxy is populated with only radial orbits. Not all density profiles can sustain such a model: Richstone & Tremaine (1984) demonstrated 
that the density must diverge at least as r"^ at small radii in order to be able to sustain a purely radial orbital structure. Models with constant 
density cores or density cusps shallower than will hence have a minimal anisotropy radius ra_miii below which Osipkov-Merritt models are 
negative at some point in phase space, and a maximum anisotropy value /J^ax above which the constant anisotropy models become negative 
at some point in phase space. For example, selfconsistent Plummer models have /3,„ax ~ and ra.min ~ f (Merritt 1985), and selfconsistent 
Hernquist models have /3niax = \ and ra.min ~ 0.202 (Baes & Dejonghe 2002). The central slope of the density J4j for the models we consider 
in this paper is steep enough such that even models with a purely radial orbital structure are positive over the entire phase space. 
For models with a completely radial orbital structure, the augmented density can be written as 

/5(*'0 = ^^, (74) 
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where the function Pr(5') is determined by the condition <63> . The distribution function is of course a degenerate function of angular 
momentum (only orbits with L = are populated), and can be calculated via the Eddington-like formula 



5(L^) p dpR d^ 



(75) 



d* - *' 

Once again, a closed expression can be obtained for the family of potential-density pairs we have discussed in this paper by transforming the 
integral to an integration with respect to the variable a. One obtains after some algebra 



327r3 



(1 - m) (1 - 4^ - 2m^ - ISg ^ 17^g + 4g^) 

+ (2-^ + £:)v/2 ^-^^t^nii 



V(2-M + £)£: 

1 + f 



(76) 



This distribution function is everywhere positive, for all values of the black hole mass /i. The tangential velocity dispersions for these models 
are of course identically zero, whereas the radial velocity dispersion can be found through equation \ifl\ . 



(^rir) = (1 - 



1 + 2r - 4r(l + r)ln 



r / 3 



1 - 4r - 8r^ + 8r^^^ 



(77) 



6 CONCLUSION 

We have described a one-parameter family of spherical models for elliptical galaxies and bulges consisting of a stellar component and a 
central black hole. All models in this family share the same stellar density profile, which has a steep cusp with a slope of — | in the centre. 
The gravitational potential for the models is a linear combination of the potential generated selfconsistently by the stars and the potential 
of a central supermassive black hole. The parameter /i, representing the mass of the black hole relative to the total mass of the galaxy, can 
assume any value between and 1 . The family therefore contains models ranging from a selfconsistent model without black hole to a model 
where the potential is entirely dominated by the black hole. We have done an extensive study of the internal and projected kinematics for 
this family of models. With the assumption of isotropy, we have calculated the intrinsic velocity dispersions, the distribution function, the 
differential energy distribution, the moments of the distribution function, the surface brightness and the projected velocity dispersions. All of 
these quantities have been expressed completely analytically, for all values of /j.. 

We have also described some extensions of the models to anisotropic orbital structures. In particular, we have considered models with 
a constant anisotropy, models with distribution functions of the Osipkov-Merritt type and models with a completely radial orbital structure. 
Also for these families of models, the distribution function and its moments can be calculated completely analytically for all values of the 
central black hole mass /x. 

We are well aware that the stellar component we have considered does not completely represent the detailed structure of the centre of 
galaxies. Firstly, it has a fairly steep density cusp, which is at the edge of the observed range in real elliptical galaxies (Lauer et al. 1995; 
Gebhardt et al. 1996). With such a steep cusp, much of the stellar mass is very centrally concentrated, resulting in an infinitely deep stellar 
potential well. The addition of a black hole potential to this stellar potential does not drastically change the global mass distribution, such that 
the effects of a black hole on the kinematical properties are probably quite conservative. For example, the distribution function is substantially, 
but not very drastically affected by the presence of a black hole, whereas this effect is much stronger for 7-models with a less steep density 
cusp (see figure 6 in Tremaine et al. 1994). The reason why we focused on the density profile J4}, and not on one of the other members of the 
family of 7-models with a shallower central density cusp for example, is that the 7 = § model is unique in the way that it allows a relatively 
simple expression of the augmented density in the presence of a central black hole. It is the only member of the family of 7-models where 
the distribution function, the differential energy distribution and the moments can completely be expressed in terms of elementary functions. 
We are currently undertaking a more general study of the properties of the 7-models with a central black hole, using both analytical and 
numerical means (Baes, Buyle & Dejonghe 2004). 

Secondly, the models presented here are spherically symmetric and isotropic, whereas few elliptical galaxies are thought to be perfectly 
spherical. Observational studies suggest that a substantial fraction of elliptical galaxies are at least moderately triaxial (Franx, van Gorkum 
& de Zeeuw 1991; Tremblay & Merritt 1995; Bak & Statler 2000). The central regions of galaxies with a central black hole are generally 
believed to be roughly axisymmetric, however. Indeed, supermassive black holes drive the shape of galaxies in a triaxial haloes toward ax- 
isymmetry by stochastic diffusion, either globally (Gerhard & Binney 1985; Norman, May & van Albada 1985; Merritt & Quinlan 1997; 
Wachlin & Ferratt-Mello 1998; Valluri & Merritt 1998), or at least the central regions (HoUey-Bockelmann et al. 2001, 2002). Still, axisym- 
metric systems have a much larger freedom in orbital structure than spherical systems. Using detailed axisymmetric dynamical modelling, 
Gebhardt et al. (2003) found that the central regions of the massive early-type galaxies generally have a significant tangential anisotropy, 
whereas less massive galaxies have a range of anisotropics. 

In spite of these critical notes, the family of dynamical models presented in this paper has the huge advantage over more complicated 
numerical models that all kinematical properties can be calculated completely in terms of elementary functions, for any value of the parameter 
IJ,. This is not only the case for rather simple kinematical properties such as the velocity dispersions, but also for more complicated kinematical 
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properties which depend in a strongly non-linear way on the potential. In particular, the distribution function and the differential energy 
distribution can be written in a compact form and in terms of elementary fimctions only. As a result, this family of models is useful for a large 

set of applications. In particular, it can be used to easily generate the initial conditions for N-body or hydrodynamical simulations, which are 
needed to investigate how black holes interact with the stellar, gaseous and dark matter components of their host galaxies. 
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APPENDIX A: EVALUATION OF THE INCOMPLETE ELLIPTIC INTEGRALS 

The expressions <51> . <56> and <59> for the surface brightness, the cumulative surface brightness and the projected velocity dispersion contain 
incomplete elliptic integrals of the first and second kind. When the modulus of these functions (the second parameter) is less than unity, 
these functions are easily evaluated numerically, and simple and efficient routines are widely available (e.g. Press et al. 2001). Sophisticated 
closed-box mathematical packages such as Maple or Mathematica easily calculate incomplete elliptic integrals for all (complex) values of the 
modulus, but most basic implementations break down when the modulus is greater than unity. For the present case, this happens when R < 1. 
A well-behaved solution for i? < 1 can be obtained with the transformtion formulae for incomplete elliptic integrals (e.g. Abramowitz & 
Stegun 1972). The relevant formulae to transform the expressions for the surface brightness profile etc. for i? < 1 are 



1" / 2 \ I R+1 ^1 _ /I 




W = E arctan Ji, + ^£=^ F arctan Jl, . (A2) 

R+lJ \ R+1 y \l R'\l 2 J ^2{R+1) \ \ R \ 2 J 

Around R — 1 either of these forms are hard to evaluate numerically, because the modulus of the elliptic integrals then approaches unity. 
Instead, it is better to use the Taylor expansion 



where 

^ arctanh f-^^ = -^\n( ^^^ ^ 0.62323. (A5) 



(A6) 
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For example, for the surface brightness profile, we obtain 
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